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It is shown that under certain conditions, three-component superconductors (and in particular 
three-band systems) ahow stable topological defects different from vortices. We demonstrate the 
existence of these excitations, characterized by a CP^ topological invariant, in models for three- 
component superconductors with broken time reversal symmetry. We term these topological defects 
"chiral GL^^^ skyrmions", where "chiral" refers to the fact that due to broken time reversal sym- 
metry, these defects come in inequivalent left- and right-handed versions. In certain cases these 
objects are energetically cheaper than vortices and should be induced by an applied magnetic field. 
In other situations these skyrmions are metastable states, which can be produced by a quench. Ob- 
servation of these defects can signal broken time reversal symmetry in three-band superconductors 
or in Josephson-coupled bilayers of s± and s-wave superconductors. 

PACS numbers: 74.70.Xa 74.20.Mn 74.20. Rp 



I. INTRODUCTION 

Experiments on the recently discovered iron pnictide 
superconductors suggest the existence of positive co- 
efficient of Josephson coupling between superconduct- 
ing components in two bands {s± state) and possibly 
more than two superconducting bands. ^ Under these cir- 
cumstances, new physics can appear. That is, frus- 
tration of competing interband Josephson couplings in 
three-component superconductors, can lead to sponta- 
neously Broken Time Reversal Symmetry (BTRS)^'^ (an- 
other scenario for BTRS states in pnictides was dis- 
cussed in Refs. 4 and 5). There, the ground state ex- 
plicitly breaks the discrete U(l) x ^2 symmetry. ^'^ Re- 
lated multicomponent states were also recently discussed, 
in connection with other materials.^ If superconductiv- 
ity in iron pnictides is described by just a two-band s± 
models, BTRS states can nonetheless be obtained in a 
Josephson-coupled bilayer of s± superconductor and or- 
dinary 5- wave material.^ Such bilayer systems can be ef- 
fectively described by a three-component model where 
the third component is coupled through a "real-space" 
inter-layered Josephson coupling. 

Due to a number of unconventional phenomena, which 
are not possible in two-band superconductors, the pos- 
sible experimental realization of three component su- 
perconductors (either with or without BTRS) recently 
started to attract substantial interest. ^'^''''^^^ These phe- 
nomena include: exotic collective modes which are dif- 
ferent from the Leggett's mode;^'^^'^^ the existence of 
a large disparity in coherence lengths even when inter- 
component Josephson coupling is very strong, leading 
to type- 1.5 regimes^ (where some coherence lengths are 
smaller and some are larger than the magnetic field pen- 
etration length^^); the possibility of flux-carrying topo- 
logical solitons different from Abrikosov vortices.^ 

This paper is a follow-up to Ref. 6 where we introduced 
new fiux-carrying topological solitons. Here we study 



in detail, these topological solitons which we term chi- 
ral GL^^^ skyrmions (chiral skyrmions for short). They 
are magnetic flux-carrying excitations characterized by 
a (DP^ topological invariant, (by contrast this invariant 
is trivial for ordinary vortices). The topological prop- 
erties, motivating the denomination skyrmion are rigor- 
ously discussed. As the terminology suggests, the soliton 
itself has a given chiral state of the Broken Time Rever- 
sal Symmetry. More precisely, different arrangements of 
the fractional vortices constituting a skyrmion carrying 
integer flux define different chirality of the skyrmion. Fi- 
nally GL^^^ refers to the physical context of the three- 
component Ginzburg-Landau theory. The thermody- 
namic and energetic (met a) stability of chiral skyrmions 
are discussed, as well as their perturbative stability. In 
scanning SQUID, scanning Hall or magnetic force mi- 
croscopy experiments, chiral GL^^^ skyrmions can (un- 
der certain conditions) be distinguished from vortices by 
their very exotic magnetic field profile. Fig. 1 shows ex- 
amples of such exotic magnetic field signatures of chiral 
skyrmions in three band superconductors with various 
parameters of the model. 

The paper is organized as follows. In Sec. II we in- 
troduce a Ginzburg-Landau model for three-component 
superconductors where phase frustration due to compet- 
ing Josephson interactions leads to Broken Time Rever- 
sal Symmetry states. The structure of the domain walls 
which are possible due to this new spontaneously broken 
Wj2 symmetry is discussed in Sec. II A. The essential con- 
cepts of the topological excitations in multi-band super- 
conductors are discussed in Sec. II C. After that, the new 
kind of topological excitations, chiral GL^^^ skyrmions, 
are discussed Sec. II D. The physical properties: (i) en- 
ergy of formation of a skyrmion versus vortex lattice, (ii) 
thermodynamical stability of the chiral skyrmions and 
(iii) their perturbative stability are investigated Sec. HI. 
In the next part. Sec. IV, the very rich interactions be- 
tween the chiral skyrmions and between skyrmions and 




Figure 1. (Color online) - Example of unusual observable magnetic field configuration of chiral skyrmions. 



vortices are investigated. The model has many inter- 
esting mathematical aspects as well. Sec. V is devoted 
to the most formal aspects and rigorous justifications of 
the physics and mathematical properties of the three- 
component Ginzburg-Landau model and the skyrmionic 
excitations therein. This section aims at a more mathe- 
matical audience. Thus, readers less interested in formal 
justification of the physics can skip these discussions, and 
go straight after Sec. IV to our conclusions in Sec. VI. 
There we conclude this paper by addressing, in more de- 
tail, the possible experimental signatures of our chiral 
GL^^^ skyrmions. 



II. THE MODEL 

In this paper we consider various realizations of three- 
component superconductivity described by the following 
three-component Ginzburg-Landau (GL) model: 



+ ^ lah\^a\^\^h? - r]ah\^a\\^h\ COs{(pb - (fa) • (2.1) 
a,b>a 

Here D = V + ieA and tjja = \^|^a\e'^'^'^ are complex fields 
representing the superconducting components. The com- 
ponent indices a, 6 take the values 1,2,3. In the partic- 
ular case of a three-band superconductor, different su- 
perconducting components arise due to Cooper pairing 
in three different bands. The bands are coupled by their 
interaction with the vector potential A and also through 
potential interactions. The coefficients rjab are the inter- 
component Josephson couplings. We also consider the 
more general case which includes bi-quadratic density 
interactions with the couplings jab- Here, the London 
magnetic field penetration length is parametrized by the 
gauge coupling constant e. Functional variation of the 
free energy (2.1) with respect to the fields gives Ginzburg- 
Landau equations 

dV 

DD^a = 2— , d, {d,Aj - djA,) = J, . (2.2) 

where V is the collection of all non-gradient terms and 
the supercurrent is defined as 



In multiband superconductors, a Ginzburg-Landau ex- 
pansion of this kind can in certain cases be formally jus- 
tified microscopically (see e.g. corresponding discussion 
in two-band case^*"). In what follows, different physical 
realizations of the model (2.1) with different broken sym- 
metries are considered. Note that in some of the phys- 
ical realizations of multicomponent GL models, some of 
the couplings are forbidden (for example on symmetry 
grounds). This can occur for intercomponent Joseph- 
son couplings, in some realizations.^^ More terms, con- 
sistent with symmetries, can be included to extend the 
GL functional. Alternatively a microscopic approach can 
provide a more quantitatively accurate picture at lower 
temperatures. However, the properties of the topologi- 
cal objects which are discussed, should then differ only 
quantitatively and not qualitatively in the framework of 
e.g. microscopic approach for a system with a given sym- 
metry (some examples how phenomenological multiband 
GL models give good results even at low temperature can 
be found in Ref. 17). 

The field configurations considered in the following are 
two-dimensional, as well as three dimensional systems 
with translation invariance along the third axis. 



Broken Time Reversal Symmetry, the U{1) 
states 



X ^2 



(2.3) 



a=l,2,3 



o=l,2,3 



For a given parameter set {aa^ Pa^Vab^ lab)^ the ground 
state is the field configuration which minimizes the po- 
tential energy. The corresponding values of |V^a|'s and 
iPa^i together with the gauge coupling e determine the 
physical length scales of the theory. The particularly in- 
teresting property of the model (2.1), is that the ground 
state can be qualitatively different from its two band 
counterparts. While in two bands systems with Joseph- 
son interactions the phase-locking is trivial (either or 
tt), the phase- locking in three bands can be much more 
involved. Indeed, competition between different phase- 
locking terms possibly leads to phase frustration. When 
> the corresponding Josephson term is minimal 
for zero phase difference, while if rjab < it is minimal 
for ipab = ^6 — V^a = TT. Now if the signs of rjab^ are all 
positive (we denote it as [+ + +]), the ground state has 

^1 = ^2 = ^3- Similarly for [H ] couplings, the phase 

locking pattern ^1 = ^2 = ^ tt. However for [+ H — ] 

or [ ], the phase locking terms are frustrated. That 

is: all three Josephson terms cannot simultaneously at- 
tain their minimal values. As a result ground state phase 
differences are neither nor tt. For example, consider the 



case = — 1, /3a = I and rjab = — 1- Symmetry under 
global U(l) phase rotations allows to set (^i = without 
loss of generality (for the below considerations). There, 
two ground states are possible (p2 = 27r/3, cps = — 27r/3 
or (p2 = — 27r/3, cps = 27r/3. The two ground states are 
each other's complex conjugate. The actual values of the 
ground state phases depend on the potential parameters. 

Note that the free energy is invariant under complex 
conjugation, (V^i, V^2, V^s) ^ i'^i^'^h'^s)^ which takes it 
to a state with different phase locking. Thus the the- 
ory has a spontaneously broken discrete (Z2) symmetry, 
called Time Reversal Symmetry. That is, the free en- 
ergy is still invariant under complex conjugation, but the 
ground state is not. By 'picking' one of the two inequiv- 
alent phase-locking patterns, the ground state explicitly 
breaks the discrete symmetry. Such states are termed 
Broken Time Reversal Symmetry (BTRS) states. 



B. Domain walls in BTRS states 

BTRS systems have topological excitations related to 
the broken discrete symmetry in the form of domain 
walls. The domain walls interpolate between domains 
of inequivalent ground states. In other words they are 
walls separating regions of different phase locking. It is 
instructive to display more quantitatively the structure 
of the ground state (or "vacuum") manifold, see Fig. 2. 
There, the potential energy is minimized with respect 
to the densities for uniform fixed phase difference 

configurations. This provides a map of the ground state 
manifold. It appears clearly that there are disconnected 
inequivalent ground states (the red and green dots). In- 
terestingly, there is not a unique path to connect inequiv- 
alent ground states with inequivalent phase locking, but 
four. The four corresponding domain walls will have dif- 
ferent line tension (energy per unit length). Note that, in- 
vestigating the vacuum manifold with fixed ground state 
densities IV^^I (at their true ground state value) provides 
a qualitatively similar picture. Namely, this approxima- 
tion preserves the positions of the minima. However, the 
actual values of J^^^^ are obviously different if | Veal's are 
held constant to the ground state, so this approximation 
does not allow one to calculate the energy of the do- 
main walls. In particular the sharp angles appear there 
for strong Josephson couplings, when the ground state 
densities are not fixed. This property is absent when 
densities are held to their actual ground state values. 



C. Flux-carrying topological defects in three 
component Ginzburg— Landau model 

As previously stated, three component Ginzburg- 
Landau model can exhibit BTRS and domain wall exci- 
tations associated with the broken symmetry. There 
are also different topological defects, associated with the 
other broken symmetries. 
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Figure 2. (Color online) - Representation of the vacuum 
submanifold (Top), for (aa,/^a) = ("1? 1) cind r]ab = —3. The 
image shows the potential energy as a function of the phase 
differences: (p2 and cps, minimized with respect to all moduli 
degrees of freedom (while cpi is set to zero by U(l) invari- 
ance associated with simultaneous change of all phases). Red 
and green dots show inequivalent ground states. And the 
lines connecting them represent four different kinds of do- 
main wall trajectory over the field manifold. Black dots are 
ground states located farther than 27v in the phase differences. 
The second line, gives a schematic representation of various 
^2 domain walls in three-band superconductors with different 
frustrations of phase angles, shown by arrows of different col- 
ors. The pink line schematically shows the phase difference 
between red and green arrow, interpolating between the two 
inequivalent ground states. 



Our main interest, here, is three- component skyrmionic 
solutions of the Ginzburg-Landau model. Here 
skyrmions are topological defects characterized by a 
topological invariant which classifies the maps M? 
CP^. In contrast to the topological invariant character- 
izing vortices {i.e. the winding number which is defined 
as a line integral over a closed path), the topological in- 
dex associated with skyrmionic excitations is given as an 
integral over xy-plane : 

(2.4) 

with ^■f' = (V^^,V^|,?/^|). A detailed derivation of this 
formula is given in Sec. V. If we have an axially symmetry 
vortex with a core where all superconducting condensates 
simultaneously vanish, then Q = 0. On the other hand, 
if singularities happen at different locations, then Q 7^ 
and the quantization condition Q = B /^o = N holds 
( ^0 being the flux quantum and N the number of flux 
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quanta). This is rigorously discussed in Sec. VB. 



Fractional vortices 

In order to understand the physical properties of the 
later introduced chiral skyrmions, it is good to remind 
oneself of the basic features of multi-component super- 
conductors and their topological excitations. The ele- 
mentary vortex excitations in this system are fractional 
vortices. They are defined as field configurations with 
a 27r phase winding only in one phase {e.g. cpi has 
Aipi = § Vipi = 27r winding while A(p2 = = 0)- To 
better illustrate their physical properties, the Ginzburg- 
Landau free energy (2.1) can be rewritten as 

^=i(VxA)^ + ^ (2.5a) 
+ E ^(Vl^al)' + + y IV'al^ (2.5b) 

a 

+ 2; T^llAPlV-lP, (2^5d) 

a,b>a 

where (fab = — are the phase differences and 
= IV^aP- The indices a, 6 again denote the differ- 
ent superconducting condensates and take value 1,2,3. 
The identity 

n n 

Yl IV^-I V6|'V(^a (V(/9a - V(/95) 
a=l 6=1 

n n 

= E E \^a\^\^b\H^Va-V^,f , (2.6) 

a=l 6=a+l 

is used to derive this expression. Here, the supercurrent 
(2.3) reads, more explicitly 

J/e = ep2A + El^a|'V^a. (2.7) 

a 

Consider now a vortex for which the phase of only one 
component changes by 27r: f S/cpa = 27r. Such a configu- 
ration carries a fraction of fiux quantum^^ 

$„=/Ad£=Mi/v^. = M$o, (2.8) 

where denotes the ground state density of i/ja^ cr is 
a closed curve around the vortex core, and ^0 = 27r/e 
is the fiux quantum. For vanishing Josephson interac- 
tions, the symmetry is [t/(l)]^ and each fractional vortex 
has logarithmically diverging energy. This can be seen 
easily in the London limit by setting tpa = const every- 
where except a sharp cutoff in the vortex core. There the 



terms (2.5d) and (2.5b) give trivial contribution to the 
free energy, so that the relevant parts now reads 

^..„„ = i(VxA)^ + ^ 

+ E^^;;^((V^-) -^-s^..j. (2.9) 

In a [[/(l)]^ symmetric model, one fractional vortex 
gives logarithmically divergent contribution to the energy 
through the term 

(2.10) 

Tc being a sharp cut-off corresponding to the core size of 
a vortex. However a bound state of three such vortices 
(where each phase a = 1, 2, 3 had 27r phase winding) has 
finite energy. Indeed such a bound state has no wind- 
ing in the phase differences. This finite-energy bound 
state is a "composite" vortex having one core singular- 
ity where + |?/^2| + iV^sl = 0. Around this core all 
three phases have similar winding A(fa = 27r. A vor- 
tex carrying one quantum $0 of flux is thus a logarith- 
mically bound state of fractional vortices. For non-zero 
Josephson coupling, fractional vortices interact linearly, 
so they are bound much more strongly. It can be seen 
that, for non-zero Josephson coupling, the phase differ- 
ence sector (2.5c) or the second line in (2.9) is a sine- 
Gordon model. There, a given fractional vortex excites 
two Josephson strings (one per phase difference sector). 
Crossections of a string, at a large distance from a vor- 
tex are sine- Gordon kinks. Such a Josephson string, has 
an energy proportional to its length. Thus for non-zero 
Josephson coupling one fractional vortex has linearly di- 
verging energy (see App. A for a detailed derivation). 
Note that the Josephson strings are different topological 
excitations than the domain walls previously discussed. 
Having linearly diverging energy, fractional vortices in- 
teract linearly. As a result an (composite) integer fiux 
vortex can be seen as a strongly bound state of three co- 
centered fractional vortices. This binding is thus much 
stronger for non-zero Josephson couplings. Because of 
their diverging energies, the fractional vortices are not 
thermodynamically stable in bulk samples: A group of 
three different fractional vortices is energetically unsta- 
ble with respect to collapse into an integer fiux composite 
vortex. Note however that under certain conditions, in a 
finite sample, they can be thermodynamically stable near 
boundaries'^ with strings terminating on a boundary. 

Note that in a London limit, magnetic field of frac- 
tional vortices is exponentially localized. However in a 
[[/(l)]^ Ginzburg-Landau model, the magnetic field of a 
fractional vortices is in a general localized only according 
to a power law and moreover can invert direction. 
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Figure 3. (Color online) - A single charge chiral skyrmion, for 3 mirror passive bands (tta, Pa) = (1, 1) and Josephson coupling 
constants r/ab = —3. Here jab = 0.8 and the gauge coupling constant is e = 0.6. Displayed quantities are the magnetic flux 
(A) and the sine of phase differences sin((/?i2) (B), sin((^i3) (C). Condensate densities (D), \ip2\^ (E) and iV^fl, (F) are 

represented on the second line. The corresponding supercurrent densities | Ji|, (G), | J2I, (H) and | J3I, (I) are displayed on the 
third line. To avoid redundant informations, the total energy density is not displayed. It qualitatively follows the magnetic 
flux shown in panel (A). 



D. Chiral three component Ginzburg— Landau 
skyrmions 

Domain- walls such as those discussed in Sec. II A can 
form dynamically in physical systems by a quench. Be- 
cause of its line tension, a closed domain wall collapses 
to zero size. From the term (2.5c), in the rewritten 
Ginzburg-Landau functional, it is clear that in order to 
decrease the energy cost associated with a gradient in the 
relative phase (pab^ the densities of the components l^/^^l, 
IV^^I should be suppressed on the domain wall. Further- 
more, on a domain wall, the cosines of phase differences 
cos{(fi) — if a) are energetically unfavorable. Indeed, by 
definition, it is where they are the farthest from their 
ground state values. As a result, if an integer composite 
vortex is placed on the domain wall, the Josephson terms 
should tend to split it into fractional flux vortices, allow- 



ing it to attain more favorable phase difference values in 
between the split fractional vortices. As a consequence of 
these circumstances, the domain wall can trap vortices. 
Recall that away from domain walls, fractional vortices 
are linearly confined by Josephson terms. 

When the magnetic field penetration length is suffi- 
ciently large (e small enough), the repulsion between the 
fractional vortices confined on the domain wall can be- 
come strong enough to overcome the domain wall's ten- 
sion. It thus results in a formation of a topological soliton 
made up of 3A^ fractional vortices, stabilized by compet- 
ing forces. Such 'composite' topological solitons are thus 
made of a closed domain wall along which there are N 
singularities in each condensate Around each singu- 
larity the phase changes by 27r. The total phase wind- 
ing around the soliton is then §V(fidl = f\/(p2d£ = 
f \/(psdi = 27tN. Therefore it carries N flux quanta. 
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Figure 4. (Color online) - A Skyrmion with Q = 6 topological charge (which implies that it carries six flux quanta and 
consists of 18 fractional vortices). Displayed quantities are the magnetic flux (A) and the sine of phase differences sin((/?i2) (B) 
sin((/?i3) (C). Condensate densities iV^il, (D), IV^II, (E) and iV^ll, (F) are represented on the second line. The corresponding 
supercurrent densities |Ji|, (G), | J2I, (H) and IJ3I, (I) are displayed on the third line. Parameters are the same as in Fig. 3. 



The (DP^ topological invariant (2.4) computed for such 
objects is found to be integer, whereas it is zero for or- 
dinary composite vortices. As a result, the composite 
configuration made out of a domain wall between two 
^2 domains stabilized by repulsion between trapped vor- 
tices, is in fact a distinct topological defect: Chiral GL^^^ 
skyrmion (chiral skyrmion for short). 

It was previously demonstrated that these topological 
defects exist and are indeed at least metastable.^ Here we 
further investigate these objects. To investigate the exis- 
tence and stability of the so-called chiral skyrmions, we 
use an energy minimization approach, using non-linear 
conjugate gradient algorithm. More details about the 
employed numerical schemes are provided in App. B. The 
topological charge (2.4) was computed numerically for all 
configurations and was found to be integer within small 
numerical errors, less than 0.1%, thus providing an esti- 
mate of the accuracy of our solutions. 

Fig. 3 shows a Q = 1 chiral skyrmion in a supercon- 



ductor with three passive bands {i.e. the quadratic terms 
have positive prefactors aa). The fact that the bands are 
passive is not important for the soliton's existence. It 
consists of three fractional vortices, each one carrying a 
fraction of magnetic fiux which adds up to a flux 

quantum ^q. Since the fractional vortices are located 
quite close to each other they cannot be distinguished 
in the magnetic field profile in this case. Single charge 
skyrmions are more difficult to obtain than higher-charge 
skyrmions in this model. As will be explained later, in- 
creasing the number of flux quanta usually makes 
the solution more stable (which contrasts with vortices 
where, in the type-II regime only N = 1 vortices are sta- 
ble). The bi-quadratic density interactions in the model 
(2.1) help to stabilize Q = 1 solutions. Single charge 
solitons are thus usually supported by bi-quadratic den- 
sity interactions. Clearly, from the density plots (panels 
(D-F)) in Fig. 3, each component has a non-overlapping 
zero (the blue spots). A feature which can be observed 
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Figure 5. (Color online) - A Q = 2 quantum soliton in a system with two identical passive bands (tta, Pa) = (1, 1) (a = 1, 2) 
coupled to a third active band with substantial disparity in the ground state densities (as, Ps) = (—2.75, 1). Josephson coupling 
constants are 7712 = 7713 = '^723 = —3. The system is in a strongly type-II regime e = 0.08, the solutions here are stable even 
in the absence of bi-quadratic density interaction i.e. jab — 0. Displayed quantities are the magnetic flux (A) and the sine 
of phase differences sin((^i2) (B) sin((/?i3) (C). Condensate densities (D), IV^II, (E) and (F) are represented on the 

second line. The corresponding supercurrent densities |Ji|, (G), IJ2I, (H) and IJ3I, (I) are displayed on the third line. 



in this regime is the strong density overshoot opposite to 
the cores (the red spots). 

Higher charge skyrmions are easily formed in many 
cases even when there is no bi-quadratic density interac- 
tion. There, the stability of the skyrmion against collapse 
of the domain wall is supported only by the electromag- 
netic repulsion and Josephson interactions. In different 
numerical simulations we quite easily constructed thou- 
sands of different skyrmionic configurations, for very dif- 
ferent parameter sets. A sample of the various skyrmions 
is given in the Figures 3-7. More regimes are given in 
the appendix App. C. For all such configurations the 
(DP^ topological charge (2.4) is integer with very good 
accuracy ( \Q/N - 1| < IQ-^ ). 

One key feature, in the Figures 3-7, is seen in the 
phase differences on panels (B) and (C). In each of these 
various regimes, the phase locking pattern 'inside' the 



skyrmion is different from 'outside', thus corresponding 
to either of the two inequivalent ground states. As 
a result the chiral skyrmions (in contrast to non-chiral) 
feature a domain wall separating the regions of different 
BTRS states. As discussed below Sec. IV B, the choice of 
one of the ground states inside the skyrmion dictates 
a clockwise versus counter-clockwise arrangement of frac- 
tional vortices, thus motivating the terminology "chiral" 
for these topological defects. 

Chiral skyrmions exhibit very unusual signatures of the 
magnetic field which can be seen from the panel (A) in 
all of the Figures 3-7 or in Fig. 1. If the bands have 
similar density, each fractional vortex carries a similar 
fraction of flux quantum. As a result, the magnetic flux 
is almost uniformly spread along the domain wall, as in 
Fig. 4. On the other hand, when the condensates have 
quite different densities, the magnetic fiux is carried non- 
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Figure 6. (Color online) - A Q = 5 quantum soliton in a system with two identical passive bands as in Fig. 5 coupled to a 
third active band with disparity in the ground state densities (asj/^a) = (—1-5, 1). Josephson coupling constants are 7723 = —3 
and ?7i2 = 7713 = 1. Here e = 0.2 and there is no density-density interaction term 'jab = 0. The system is shaped as a pentagon 
deformed by the vortices of the strong active band carrying larger fractions of flux quantum. Displayed quantities are the same 
as in the previous pictures, e.g. Fig. 5. 



uniformly by fractional vortices in different condensates. 
Consequently, the magnetic flux is inhomogeneously dis- 
tributed along the soliton. This can be seen in Fig. 5 
where the third component carries a great fraction of the 
flux. The remaining fraction of flux is spread along the 
components having less density. The overall configura- 
tion can easily be mistaken for a vortex pair in such a 
superconductor. For higher topological charge, the same 
system exhibits geometric structures (a pentagon as in 
Fig. 6) where the vertices are occupied by the fractional 
vortices of the band with bigger density. There again, 
geometrical arrangement of apparent vortices is a very 
typical signature of the chiral skyrmions. 

Among possible observable signatures of chiral 
skyrmions, is the varying fraction of magnetic flux car- 
ried by fractional vortices, as in Fig. 7. There, the mag- 
netic field exhibits spots of different magnitude, larger 
spots associated to the two similar bands with more den- 



sity while the small spots are associated with the active 
band. 



E. Chiral multi-skyrmions 

Besides having non trivial (DP^ topological invari- 
ant (2.4), the chiral skyrmions in three component 
Ginzburg-Landau theory with BTRS have a given chiral- 
ity. Namely, there is a difference whether one or the other 
broken state is 'inside'. Here we report bound states 
of chiral skyrmions with opposite chirality which can be 
called multi-skyrmions. More precisely a bound state of 
a skyrmion with a given chirality, carrying some topolog- 
ical charge say Qi and a skyrmion with the opposite chi- 
rality carrying Q2, see Fig. 8. There the inner skyrmion 
has a smaller charge than the outer one, Qi < Q2 since 
the chiral skyrmion's size is controlled by the number of 




Figure 7. (Color online) - A Q = 5 quantum soliton in a system with within the same parameter set as in Fig. 6 apart from 
(asj/^s) = (—0.5, 1). Displayed quantities are the same as in the previous pictures, e.g. Fig. 5. 



enclosed quanta. The bigger is the difference between 
Qi and Q2, the weaker is the interaction between the 
two chiral skyrmions. Conversely, as Qi Q2 the chi- 
ral skyrmions interact progressively more strongly. For 
very close values of Qi and Q2 the chiral skyrmions falls 
into each other's attractive basins and the domain walls 
annihilate. This allows decay to ordinary vortices. 

Note that "opposite chirality" should not be confused 
with opposite flux, i.e. these objects have opposite chi- 
rality because they interpolate between two different 
ground states. In that respect in the BTRS case, an ad- 
ditional ^2 topological charge like those of ordinary do- 
main walls can be attributed to skyrmions. However hav- 
ing opposite ^2 topological charges does not mean that 
these objects represent a skyrmion and an anti-skyrmion. 
This is because they have similar signs of Qi and Q2 
charges as well as similar signs of the total phase wind- 
ing in the local U{1) sector. That is, they carry magnetic 
flux in the same direction. For a given skyrmion one can 
construct an anti-skyrmion from similar number of anti- 
vortices. Using anti- vortices changes the overall phase 



winding and thus the direction of carried flux. As will be 
clear from the discussion below, an anti-Skyrmion with 
the same Z2 charge as a Skyrmion will also have frac- 
tional vortices arranged in a different order. 



Similarly, there exist also "Russian nesting doll"- 
like multi-skyrmions made of larger number of alternat- 
ing skyrmions of opposite chiralities. Such a multiple 
skyrmion can be seen in Fig. 9 which shows tri-iing solu- 
tions of skyrmion with alternating chiralities. This kind 
of numerical solution is quite easily obtained given a good 
initial guess. However this configuration can also sponta- 
neously form from 'collisional dynamics' of energy mini- 
mization of an initial configuration of closely spaced or- 
dinary vortices. This indicates that formation of multi- 
skyrmion solutions does not in general require fine tun- 
ing. Instead these solutions have a substantial "attrac- 
tive basin" in the GL energy landscape indicating they 
could also be observed in three component superconduc- 
tors with Broken Time Reversal Symmetry. 
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Figure 8. (Color online) - A Q = 11 quantum mi^/tz-soliton in a system with three identical passive bands as in Fig. 13. The 
current soliton is not made out of one but two stabilized domain walls thus being a homogeneous 6z-ring configuration. Panels 
(B) and (C) clearly display the alternating different ground- states. Since the three bands are identical, the magnetic field 
rather homogeneously spreads all along the solitons. Displayed quantities are the same as in the previous pictures, e.g. Fig. 5. 
Note that while going counterclockwise along the outter ring, the fractional vortices have order band-" 1,2,3" . For the inner 
ring they are ordered as band-" 1,3,2" . The origin of this is discussed in Sec. IV B 



III. PHYSICAL PROPERTIES OF CHIRAL 
SKYRMIONS 



It is important to know the energetic properties of 
skyrmions compared to ordinary vortices, as well as their 
stability properties. Indeed if skyrmions are thermody- 
namically stable and form as the ground states in mag- 
netic field, their experimental signatures are straightfor- 
ward to detect. However, if they form as states with 
higher energy than e.g. a vortex state, they are only 
metastable. When they are metastable states, skyrmions 
are protected against decay by an energy barrier. The 
height of this barrier depends non-trivially on the pa- 
rameters of the potential and on the number of enclosed 
flux quanta. Metastable chiral skyrmions could be pro- 
duced by quenching the system under applied magnetic 
field. In this section, we discuss these aspects. 



A. Energy of Chiral skyrmions vs vortices 

For vanishing bi-quadratic density interaction cou- 
plings {i.e. Jab = 0), in all the regimes which we in- 
vestigated, chiral skyrmions are always more expensive 
energetically than vortices. However, as suggested in 
Ref. 6, bi-quadratic density interaction decreases the en- 
ergy of chiral skyrmions relative to that of vortices. For 
sufficiently strong bi-quadratic density interaction chiral 
skyrmions are ground state excitations i.e. energetically 
cheaper than vortices and, for certain parameters, ther- 
modynamically stable. 

The energy properties of the chiral skyrmions are dis- 
played on the left panels of Figures 10-11. There, the 
energy per flux quantum of a given configuration is given 
in units of the single quantum fiux carrying ground state. 
Namely E{N)/[NE{N = 1)] is represented as a function 
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Figure 9. (Color online) - A tri-iing chiral skyrmion. The configuration carries total charge Q = 36 in a system with two 
identical passive bands = (0^2, A) = (1, 1) coupled to a third active band with (0^3, /^s) = (—1, 1). Josephson coupling 

constants are 7723 = —3 and 7712 = 7713 = 1. e = 0.7. Panels are the same as usual, e.g. Fig. 5 



of the number of flux quanta. The corresponding en- 
ergies are sublinear functions of enclosed flux quanta for 
all solutions with N > 2. This means that the energy 
cost per flux quantum decreases as N grows. 

Two different regimes can be distinguished. If a con- 
figuration has E{N)/[NE{N = !)]>! (where E{N = 1) 
is the energy of a single vortex), then it is energetically 
preferable to have N isolated type-II integer flux vor- 
tices. As discussed below, there, skyrmions should be 
understood as metastable objects. That is, they can 
decay into type-II (composite) vortices, e.g. in case of 
strong enough perturbations. On the other hand, when 
E{N)/[NE{N = 1)] < 1, then isolated vortices are no 
longer energetically preferred over a skyrmion. In the 
first case, (corresponding to the upper curves of Fig. 10), 
chiral skyrmions can exist as meta-stable excitations. In 
the second situation (the lower curves of Fig. 10), chiral 
skyrmions could form as true ground state topological 
excitations. Note also that there is a regime where lower 
charge skyrmions are more expensive than type-II inte- 



ger vortices, while higher charge ones are cheaper (see 
Fig. 10). In the regimes where there is density-density 
interaction, even the smallest skyrmions with Q = N = 1 
can be energetically cheaper than vortices. 

The relative cost of including an additional flux quan- 
tum into a chiral skyrmion is evaluated by computing 
^''^E{N=i)~'^'' • When this quantity is less than one, it is 
globally beneficial to merge an additional flux-quantum- 
carrying object with a skyrmion. It is displayed on panel 
(B) of Figures 10-11. Note that it does not tell about 
the real work the system has to provide for bringing 
the isolated single quantum defect from infinity into the 
skyrmion, but only on global cost or benefit. 



B. Thermodynamical stability of Chiral skyrmions 

The first critical field is defined as the applied mag- 
netic field at which the formation of a single flux car- 
rying defect (vortex or skyrmion) becomes energetically 
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Figure 10. (Color online) - Energies per flux quantum of the skyrmions carrying N flux quanta. The energy is given in units 
of the energy of the energetically cheapest (either vortex or skyrmion) single quantum excitation (A). Middle panel (B) shows 
^^^(A^i)""*^^ ^ function of the number of flux quanta. When this quantity is less than one, it is energetically preferred to 
have a A/'-quantum skyrmion than having a (N — l)-skyrmion plus one isolated vortex. The criterion for thermodynamical 

stability — |Fgs|, where the condensation energy is Fgs = F({ipa), 0), is shown on the right panel (C). Here, the dependence 
of the solutions on 7 and N is investigated, while the gauge coupling is fixed at e = 0.3. Other parameters are (tta, Pa) = (1, 1) 
and riab = —3. Colors and symbols associated to different values of jab (shown on the picture) are the same over three panels. 
Note that the reason why curves with high jat terminates for smaller N, is that the size of the skyrmion becomes comparable 
to the size of the numerical domain. To avoid any finite size effect, we chose to skip the corresponding points. 



favorable. It is defined in analogy with the first critical 
field for ordinary vortices Hd = E^/^di where Ed and 
are the energy and magnetic flux of the topological 
defect Ed = J {J^i^l^a, A) - J^^s) and J^gs = ^((V^a), 0) is 
the ground state energy. I.e. it is energetically preferred 
to form a topological defect carrying flux ^d in external 
field Hq if the Gibbs free energy Ed - ^dHo < 0. The 
external field Hq should be smaller than the thermody- 
namical critical mao: netic field Hd = 2^/J^{0,0) - J^gs- 
The criterion for thermodynamical stability is investi- 
gated on the right panels (C) of Figures 10-11. For all 

these regimes, — {Tcsl < 0. In all displayed cases, 
skyrmions satisfy this criterion. That means that un- 
der certain conditions they can be induced by an applied 
external field. 



C. Perturbative stability of Chiral skyrmions 

Chiral skyrmions can appear as thermodynamically 
stable ground states or metastable states in supercon- 
ductors with Broken Time Reversal Symmetry. In this 
work, they are obtained by minimizing the energy. Con- 
sequently, they are always minima (at least local) of the 
free energy landscape. When the chiral skyrmions are 
metastable states they are protected against decay into 
type-II vortices by a finite energy barrier. The analysis 
carried out in this subsection concerns the metastable 
solutions. In all the regimes which we considered, 
metastable chiral skyrmions are found to be very robust. 
They are easily formed during the energy minimization, 
e.g. in closely spaced groups of vortices. The energy bar- 
rier preventing them from decay to type-II vortices is 
typically quite high. Although difficult to quantify, it 
is interesting to have a qualitative insight into the be- 
haviour of metastable skyrmions against fluctuations. 



One possible approach to study the stability of 
skyrmions is the linear stability analysis which consists of 
applying infinitesimally small perturbation to the fields, 
and investigating the eigenvalue spectrum of the (linear) 
perturbation operator, on the background of a given so- 
lution. When the background solution is (meta) stable 
all infinitesimally small perturbations are positive modes 
and thus can only increase the energy. As a result lin- 
ear stability analysis cannot tell anything especially in- 
teresting about the properties of skyrmions. A strong 
perturbation should cause a decay of a metastable chiral 
skyrmion to ordinary vortices. Here, the stability is in- 
vestigated numerically by perturbing the chiral skyrmion 
by white noise. This allows one to investigate the full 
non-linear response where the meaningful information 
belongs. The white noise applied to all degrees of free- 
dom, is generated as follows 

V^, = V^io)+Pmax(|V^|)/i^(x,^), 

A, = Af^ +Pmax(|A|)/if(x,y). (3.11) 

Here denotes the fields of the initial skyrmionic state, 
P is a ratio giving the relative magnitude of the pertur- 
bation with respect to the maximal amplitude of a given 
field of the initial state. /i^(x,7/), and iif{x^y) are (in- 
dependent) random functions of the space. They satisfy 
l/i^l < 1 and I /if I < 1. As a result all fields initially 
receive a noise whose relative amplitude is P. The per- 
turbation has very large field gradients since it is applied 
locally on the mesh. After applying noise the system 
is then relaxed using the same minimization scheme as 
for constructing the skyrmions. Despite the strong field 
gradients, if the white noise does not exceed a certain 
threshold, the configuration relaxes back to the initial 
chiral skyrmion solution. This can be seen from the up- 
per panel of Fig. 12. The noise was gradually increased, 
confirming that indeed, a sufficiently strong perturbation 
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Figure 11. (Color online) - Energies per flux quantum of the chiral skyrmions, in the units of the energy of the energetically 
cheapest (either vortex or skyrmion) single quantum excitation (A). Curves with same color and symbols on different panels 
have same parameters. The middle panel (B) shows that it is always beneficial (within a parameter range) to have a higher 
charge skyrmion than a lower charge one plus an isolated one quantum vortex. The criterion for thermodynamical stability of 

A^-quantum solitons — |J^gs|, where J^gs = J^(('0a),O) is the condensation energy (C). The dependence of the solutions 
on e and N is investigated, for a strength of the density-density interactions jab = 0.8. Other parameters are the same as in 
Fig. 10. Here again curves are truncated when the soliton's size becomes comparable to the numerical domain. 



drives the metastable solution over the barrier, in the en- 
ergy landscape, thus leading to its decay to ordinary vor- 
tex solutions as shown on the bottom panel of Fig. 12. 
The precise value of the relative amplitude required to 
destabilize a given chiral skyrmion, obviously depends on 
the parameters of the Ginzburg-Landau functional and 
on the number of flux quanta of the solution. 

As expected, if a perturbation is strong enough, the 
metastable chiral skyrmion decays to the configuration 
with less energy, i.e. isolated type-II vortices. The ob- 
served behaviour confirms the expectations from energy 
arguments Sec. Ill A. Moreover, the deeper in the type- 
II regime, the less breakable are the skyrmions. One of 
the easiest ways for a skyrmion to decay is to deform it 
enough so that the domain wall self intersects. The con- 
figuration then can decay to skyrmionic configurations 
with lower Q which are less stable and can further decay 
into integer vortices. 



IV. INTERACTIONS OF CHIRAL SKYRMIONS 

The analysis of the energetic properties of chiral 
skyrmions suggests they should have quite non trivial 
interactions. Generally, the energy per fiux quantum de- 
creases with the topological charge (see e.g. Fig. 10). In 
some cases it is also preferable to absorb isolated vor- 
tices into a skyrmion, i.e. the energy of an A/'-quantum 
vortex is less than that of an {N — l)-quantum vortex 
and an isolated vortex. In those cases, the interaction 
at short range should be attractive. On the other hand, 
they exist in regimes where vortices usually exhibit repul- 
sive interaction (type-II or even type- 1.5). Moreover, the 
lack of axial symmetry and complicated internal struc- 
ture featuring fractional vortices can provide very non- 
trivial contribution to the interaction of skyrmions in 
BTRS superconductors. 



A. Chiral skyrmion— vortex interaction 

Chiral skyrmions can have very non trivial, non- 
monotonic interaction with vortices. As seen from the 
numerically obtained solutions shown on Fig. 10 and 
Fig. 11, in applied field, chiral skyrmions can be either 
ground states (for a given phase winding) or represent 
metastable states. For some regimes, as seen from the 
middle panels of Fig. 10 and Fig. 11, a vortex placed suf- 
ficiently close to a chiral skyrmion should be absorbed in 
the domain wall and split into fractional vortices, thus in- 
creasing the charge of the skyrmion and then decreasing 
its energy per flux quantum. Consequently, the interac- 
tion is expected to be attractive at short range. Indeed, 
as we observe in numerical calculations, if vortices are 
placed close enough to a domain wall, they are easily 
trapped to form a skyrmion of larger topological charge. 
However the long range forces between skyrmions and 
vortices can be repulsive. This is clearly seen from the 
existence of stable configurations where a number of in- 
teger flux vortices are confined within a chiral skyrmion, 
as shown on Fig. 13. That figure demonstrates that there 
is a repulsion between inner "ordinary vortices" , and the 
fractional vortices comprising the chiral skyrmion, which 
follows from (i) the stability of the configuration and 
(ii) the fact that the type-II vortices visibly stretch the 
skyrmion. Thus the interaction here is non-monotonic, 
being long range repulsive, but short range attractive. 

The repulsive long-range skyrmion-vortex interaction 
follows from the following considerations. In the ground 
state a vortex is an axially symmetric object with all 
phases winding around the same core. Thus in the type- 
II limit its energy and long-range interactions are dom- 
inated by the supercurrent J term in (2.5a). At long 
separations when linearized theory applies, the inter- 
action between a skyrmion and a vortex is dominated 
by this current-current J-mediated interaction, result- 
ing in repulsion. The attractive interaction at short dis- 
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Figure 12. (Color online) - Relaxation of a randomly perturbed chiral skyrmion. Displayed quantities are the energy density, 
Im('0i'02) and The parameters are the same as in Fig. 3, but vanishing bi-quadratic density interactions ^ah — and 

e = 0.3. Thus it is only met a- stable. The snapshots show the state of the system at different stages of the energy minimization 
algorithm after the applied perturbation. On the top panel, a Q = 6 chiral skyrmion with initial white noise of 70 % of the 
ground state values. The configurations relaxes to a chiral skyrmion. On the bottom panel, a perturbation of a metastable 
charge Q = 3 soliton with an initial noise P — 0.8. Here, the noise is strong enough to break up the domain wall. The soliton 
thus relaxes to ordinary type-II vortices (one can clearly see the disappearance of the domain wall between blue and red area 
in the middle row). The last snapshot in the lower configuration does not represent a stationary configuration: the vortices 
repel each other and are in process of drifting apart. 



tances is a nonlinear effect where split fractional vortices 
in a Skyrmion can deform a vortex by "polarizing" it. 
i.e. they can split its constituent fractional vortices thus 



inducing "dipole"-like interactions. This interaction at- 
tracts the vortex so that it merges into the skyrmion. 
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Figure 13. (Color online) - A Q = 9 quantum configura- 
tion of mixed vortices and skyrmions in a system with three 
identical passive bands as in Fig. 12. This configuration is 
made out of a skyrmion surrounding two ordinary vortices. 
It is known, from energy considerations that interaction is 
short range attractive. Interaction with vortices deforms the 
skyrmions. This shows that it is long range repulsive. 



B. Skyrmion— skyrmion interaction 



In contrast to ordinary vortices in Ginzburg-Landau 
theory, chiral skyrmions do not exhibit rotational sym- 
metry. An important consequence is that inter-soliton 
interactions should in general depend on the relative ori- 
entation of the solitons. First, note that the orientation 
and position of a soliton can be described by the posi- 
tion of the fractional vortices. The shape of a soliton, 
including the positions of the constituting fractional vor- 
tices is determined by energy minimization. The energy 
of the skyrmion is invariant under overall rotation and 
translation. 

Finally note that there are two orders in which the frac- 
tional vortices can be arranged. Going counter-clockwise 
along the domain wall, the vortices can be ordered 1, 2, 3 
or 1,3,2. We denote this order o = e^bc, e being the 
Levi-Civita symbol and a, 6, c are the band indices of the 
fractional vortices. For a skyrmion carrying integer flux, 
o = ±1 (note that this ordering closely relates to the con- 
cept of chirality). As illustrated in Fig. 14 (a), a system 
of two solitons is thus described by the distance between 
them i?, their relative orientation v together with the 
ordering (chirality) of each individual skyrmion. 
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Figure 14. (Color online) - Panel (a) shows a schematic pic- 
ture of how soliton interactions are computed. This figure 
shows the interaction between two single quanta solitons, each 
consisting of three fractional vortices shown in green, blue and 
red. This generalizes easily to larger solitons. One soliton (x) 
is placed in the origin, while the second (y) is placed at a 
distance R at an angle v. Consequently, as v is varied, the 
relative orientation of the solitons changes. Case 1 shows a 
system of two solitons with identical chiralities (same ordering 
o), while case 2 shows two solitons with opposite chiralities 
(different o), although the mirrored soliton is not necessar- 
ily stable. A schematic comparison of solitons with different 
ordering o is displayed on panel (b). In the case (2), the 
gradients in phase difference due to the fractional vortices 
naturally interpolate between two states. For this reason, 
the case (2) is energetically preferable over (1) and it was ver- 
ified numerically. Finally, panel (c) gives a schematic view of 
the merging of two single quanta solitons. In order to merge, 
they should have same ordering but opposite orientation. 



1. Chirality of skyrmions: inequivalence of left- and right- 
handed solutions 

In general, for a chiral skyrmion, the energy is not in- 
dependent of the ordering o. For a given ground state 
outside of a skyrmion, the system allows only one partic- 
ular ordering o of the fractional vortices in the skyrmion. 
The mechanism that gives rise to this behaviour is illus- 
trated in Fig. 14 (b): For a given external phase- locking 
pattern (a state), only a particular ordering o gives 
the opposite state inside. In the illustration the two 
solitons (case 1 and 2) differ in the ordering of the frac- 
tional vortices (represented by red blue and green dots 
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with band index 1,2,3 respectively) - the corresponding 
phase configurations are shown by the arrows. Thus, the 
ordering of the first one (case 1) is o = 6132 = — 1 while 
the ordering of the second (case 2) is o = 6123 = +1. 
Now for a same given ground state outside both solitons, 
the phase-locking inside is determined consistently with 
the phase gradients of each fractional vortex. In the first 
case, it results in a phase arrangement inside the soliton 
that is not a ground state. However, in the second case, 
the state obtained inside is a different Z2 ground-state. 
As a result, there is a synergy effect where the phase gra- 
dients due to the fractional vortices go from one Z2 state 
to another. Therefore o = +1 is energetically cheaper 
than o = —1 for which the inner phase locking is the 
farthest from the ground state. This is indeed confirmed 
in our numerical simulations where a skyrmion o = — 1 
decays into a skyrmion o = +1. Thus the ordering of 
the fractional vortices does matter in BTRS supercon- 
ductors. It results in the discrimination of one ordering. 
This further motivates the terminology chiral. 



2. Numerical calculations on inter- skyrmion forces 

As illustrated in Fig. 14 (a), inter-soliton forces are 
computed according to the following procedure. First, 
the structure of the soliton is determined by uncon- 
strained energy minimization, thus determining the ac- 
tual position of the fractional vortices constituting the 
skyrmion. Then two skyrmions (x and y in Fig. 14) are 
place at a distance R and a relative orientation v. There, 
the energy is minimized with respect to all degrees of 
freedom, except the position of the singularities of each 
fractional vortex. As shown in Fig. 14 (a), the energy 
is computed for every distance and relative orientation 
R and v. While allowing computation of long-range in- 
ter soliton forces, this procedure has an important lim- 
itation. It does not take into account one of the non- 
linear effects: Deformation of interacting solitons in the 
form of changes of the position of the fractional vortices. 
However, this is primarily a problem at short separation, 
where the deformation is generally the strongest. 

Fig. 15 (a) shows the interaction energy of two single 
quanta skyrmions, identical to the one in Fig. 3. From 
Fig. 11 it is clear that the energy per flux quanta de- 
creases with the number of flux quanta. For the solitons 
to merge, they need to have opposite orientation, see 
Fig. 14 (c). The computed interaction energy. Fig. 15 
(a), is indeed consistent with this picture. When the rel- 
ative orientation, v is not optimal i.e. v tt, the solitons 
exert a torque on each other, so that they attain this op- 
timal orientation. Then, an attractive channel opens in 
the potential, allowing them to get closer where nonlinear 
effects are strong, ultimately leading to a merger. 

The interaction energy of a slightly more complex soli- 
ton is shown in Fig. 15 (b). There, each skyrmion carries 
two flux quanta {i.e. their topological charge is Q = 2). 
The parameters are the one of Fig. 5, from which we know 
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Figure 15. (Color online) - Panel (a) displays the interaction 
energy of two single quantum skyrmions. One soliton is placed 
at the origin, the interaction energy is plotted as a function 
of the position and relative orientation of the second soliton. 
The interaction energy is maximal when v — while it is 
minimal for the opposite orientation, v — 71. The strength of 
the interaction decreases with the separation R. The model 
parameters are the same as in Fig. 3. Panel (b) shows the 
interaction energy of two Q = 2 quanta solitons, for the same 
parameters as in Fig. 5. Note that the skyrmion has two- fold 
symmetry (it is invariant under global rotations of tt). The 
minimum energy is found for the relative orientation v^7t/2. 



that superconducting components are not identical and 
that the skyrmion is more or less elliptic. Global orien- 
tation of the skyrmion is chosen so that when v = the 
major axis of both solitons lie along the horizontal axis. 
Note that these skyrmions are not only invariant under 
global rotation by 27r, but also by tt. Within the nu- 
merical accuracy, the inter-skyrmion interaction is always 
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repulsive. Note that this approach can accurately deter- 
mine the interaction only at sufficiently long distances. 
Indeed, by fixing the positions of the fractional vortices, it 
assumes that the skyrmions are almost-rigid bodies. The 
relative position of singularities in each fractional vortex 
is fixed once for all, but the fields can deform around this 
rigid 'skeleton'. This neglects the possibility of mutually 
induced deformations of the 'skeleton', which can open 
an attractive channel. Since our "almost-rigid body" ap- 
proximation holds only at large enough distances, short 
range data are irrelevant and not displayed in Fig. 15. 
We also derive general long-range intersoliton forces in 
the more formal framework of Sec. VC. In Sec. VE the 
formal long-range interactions are applied to the partic- 
ular case of a BTRS superconductor. The predictions 
derived there are consistent with the numerical results 
presented in this section. 



V. MATHEMATICAL ANALYSIS OF 
LONG-RANGE INTERSOLITON FORCES 

The model considered in this paper has many prop- 
erties that are interesting from a formal, mathematical 
point of view. In this section we show how, by re-writing 
the free energy in terms of gauge-invariant fields, we can 
identify a hidden topological charge, associated with the 
topology of the complex projective space CP^, and devise 
a mathematically satisfactory scheme for deducing the 
nature (attractive or repulsive) and range of the dom- 
inant force between well-separated solitons (either vor- 
tices or skyrmions). For generic parameter choices, the 
final step in this scheme (finding the spectrum of a sym- 
metric real matrix) must be done numerically, but there 
are several symmetric cases and parametric limits where 
all calculations can be completed explicitly. After treat- 
ing the general case, we consider two such special cases, 
both of potential phenomenological interest. 



A. Reduction to a supercurrent coupled CP^ ^ 
model 

In this section we consider a general k component GL 
model, with no restriction on the potential terms V. The 
k complex fields tjja may be collected into a complex k- 
vector ^ : M ^ C^, where M = denotes physical 
space. It is convenient to use polar coordinates on by 
defining 



^ =: pZ 



(5.12) 



where p = V^t^ > and ZtZ = 1. Let tt : C^\{0} ^ 
(^pk-i (denote the canonical projection which takes a 
point in to the complex line through containing 
that point, and for any X G C^, X ^ 0, denote by [X] 
its projective equivalence class (so [X] = 7r(X)). By 
gauge invariance, the potential V{^) can actually de- 
pend only on p and [Z] G CP^~^, the projective equiv- 



alence class of Z, or, equivalent ly, of ^. Let ^ = tt o ^. 
This is a CP^~ ^-valued field which maps each p G M 
to [^(p)] = [Z{p)] e CP^"^ By construction it is, like 
p, gauge invariant. We may rewrite the free energy en- 
tirely in terms of the gauge- invariant quantities p, ^ and 
J = eIm(^'''P)^), the total supercurrent. To do so, it 
is convenient to think of the gauge field A and the su- 
percurrent J as one-forms rather than vector fields (so 
we use the metric on physical space M = to "lower 
the indices" on vectors and J*). In this language, the 
covariant derivative of ^ is, likewise, a one-form 



P)^ = d^ -h ieA"^ 

with values in C^. 

On C^\{0}, let us define the real one-form 



-Im- 



XtdX 

JxF' 



(5.13) 



(5.14) 



where X = (Xi, . . . , Xj.) is a global coordinate on C^\{0} 
and dX = (dXi, . . . , dX^) are the corresponding holo- 
morphic one- forms. Then the total supercurrent is 



(5.15) 



where ^*z/ denotes the pullback of u e Q^{C^\{0}) to 
M by the map ^ : M ^ C^\{0}. In less compact 
notation, this is the one-form on M whose dx* compo- 
nent is — p~^lm'^^ di"^ . It follows that the magnetic field 
(thought of as a two-form) is 



B = dA 



d(^*i^) 



J 



(5.16) 



It is a general fact that the exterior differential opera- 
tor d commutes with pullback of differential forms, so 
d(^*i^) = ^*(di^). Note that du is a closed two- form on 
C^\{0}. Let h denote the Fubini-Study metric on CP^"^ 
with constant holomorphic sectional curvature 1, and uj 
denote its associated kahler form. Then the pullback of 
by TT : C^\{0} CP^~^ is, like dv, a closed two- form 
on C^\{0}. In fact, u is defined^^ by the requirement 
that 



7r*a; = 2du. 



(5.17) 



Hence 



d(^*z/) = ^*{diy) = -^*(7r*cj) 



i(^o^)*c^ = i$*c^, (5.18) 



and so 



e V2 e VP 



(5.19) 



Similarly, we may rewrite |P^p entirely in terms of 
the gauge invariant quantities p, <I> and J. From (5.15), 
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we see that 



pZ. (5.20) 



Let ei , 62 denote an orthonormal frame on M (for ex- 
ample Ci = and Xi = dZ{ei) G TzS'^^'^. Then 
Re(Z'''Xi) = since Xi is tangent to the unit sphere in 
at Z. Hence 

i 

+ 2Im{X} Z)p' 



ep^ 



e^p^ 



|J|^+p^^(|X. 



(5.21) 



since Im(X/Z) = v{Xi) = (^*i/)(ei). Consider 7r*/i, the 
puhback by tt of the Fubini-Study metric on CP^~^ to 
C^\{0}. Given any tangent vector X e TzS'^^~^, 

(7r*/i)(X,X) = /i(d7rX, dTrX) = cj(d7rX, idyrX) 
= cj(d7rX, dTTzX) = 7r*cj(X, iX) 
= 2du{X, iX) = 4(|Xp - i^(X)2) (5.22) 

where we have used the fact that tt : C^\{0} CP^"^ 
is holomorphic (so dyr commutes with i). Hence 

i i 

= ^^7r*/i(d^ei,d^ei) 

i 

= i^Md*ei,d$e,) = J|d$|2, (5.23) 

where |di>| denotes the norm of the hnear map d$p : 
TpM — > r$(p)CP''~-'- with respect to the metric h. Sub- 
stituting (5.23) into (5.21), one sees that 



\D^\ 



Id$l 



(5.24) 



Finahy, we obtain an expression for the total free energy 

I M 




V{p,^)). (5.25) 



The above expression for F is valid for any number 
of condensates /c, and for all field configurations where 
^~^(0) C M has measure zero, i.e. where the set of 
points in physical space at which the condensates ijja all 
simultaneously vanish is negligible. This condition holds 
for skyrmions (^~^(0) is empty), and for (multi-) vortices 
(^~^(0) is finite), so we can use (5.25) for questions in- 
volving either type of soliton, though one should note 
that, for vortices, the CP^~ ^-valued field ^ is undefined 
at the finite collection of vortex positions. 

In the special case /c = 2, we may identify CP^~^ with 
the unit two-sphere S'^ , by mapping [^1,^2] G CP^ to 
the point on S'^ with stereographic coordinate Z2/Z1, so 
that $ can be interpreted as being two-sphere valued. 
The kahler form uj coincides with the area form on S'^ 
under this identification, so that the expression for F 
(5.25) reduces to the decomposition in Ref. 22. In the 
general k case (which was previously discussed, in some- 
what different mathematical language, in context of an 
SU{N) model in Ref. 23), the field <1> takes values in 
£^pk-i^ which we cannot identify with any sphere. 



B. Flux quantization and the topological charge 

In order for a configuration on M = to have finite 
total energy, ^ and p should tend to constants ^0 ^ 
CP^~^, po G (0.00), and J should tend to as \x\ 00. 
It follows, from (5.19) and Stokes's theorem, that the 
total magnetic fiux of a finite energy configuration is 



/ P = 1 / $*cj =: 



27r 



Q(^), (5.26) 



which is a homotopy invariant of the map ^ : M ^ 
CP^~^, because u is closed. In the case /c = 2, Q is the 
winding number of the map ^ : M ^ S^. For /c > 2, Q is 
still an integer, but its geometric interpretation is more 
subtle: the image of M under ^ is homologous to Q{^) 
copies of the generator of H2{CP^~^). This gives an 
alternative interpretation of Q, to augment the physical 
interpretation, described in Sec. II C, of the magnetic fiux 
being carried by an integer number of sets of k fractional- 
fiux vortices. 

It is straightforward to give an integral formula for 
Q{^) in terms of the original condensates ^, using the 
fact that 7r*a; = 2dz^: 

<I>*cj = (tt o ^)*a; = ^*(7r*cj) = 2^*dz/ 

_ 2 /dZt AdZ ZtdZ AdZtZ\ 

-7 V \z\^ ^ ; 

d^pd^"^ Ad^ + ^M^ Ad^"^^) . (5.27) 



Hence 



Q(*) = / 



2^|^|4 



(5.28) 
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One should note that the flux-quantization condition 
(5.26) and the integral formula for the topological charge 
Q above are valid only for fleld conflgurations for which 
^ never vanishes. Note that flux is also quantized for 
ordinary vortices, for which ^ vanishes, but then it is no 
longer associated with the the topological charge Q, but 
with a topological charge associated with the total 
phase winding at spatial inflnity. This expression for Q 
can be easily discretized for use on a numerical lattice. 
Comparing Q with the total number of flux quanta gives 
a convenient way of distinguishing between vortices and 
skyrmions numerically. 



C. Long-range intersoliton forces 



Let Ei^ i = l,...,2/c — 1 be an orthonormal basis of 
eigenvectors of Hp^ with corresponding eigenvalues mf > 
0. Then we can expand (cr, F) G Tp^P relative to this 
basis 



2/c-l 



whereupon we obtain 



(5.33) 



FHn = l [ (^(|dJ|^ + eVg|J|^) 

+ ^(|dai|2+mfaf)l. (5.34) 



The key to understanding long-range forces between 
solitons is to identify the point sources which replicate, 
in the linearization of the fleld theory about the vacuum, 
the asymptotic flelds of an isolated soliton.^^ Assuming 
that the vacuum is not ^ = 0, we can use the gauge- 
invariant variables p, <l>, J, and expression (5.25) for this 
purpose. So, let the vacuum {i.e. minimum of V) occur 
at p = po, ^ = ^0- To identify the linearization of the 
theory about this vacuum, we set p = po + ^r^ ^ — *^o + ^7 
where Y G T<i>qCP^~^, and expand F to quadratic order 
in the small quantities a, Y and J: 



1, 



Idal- 



-HeSS(p,,<^,o)((cr,y'),(cr,r)) 



+ 



2eV^ 



(|dJp 



(5.29) 



where HesS(pQ^<^Q) is the Hessian of the function V 



(0, oo) X CP 



k-l 



about its minimum (po, ^^o), which 



we now deflne. Let P = (0, oo) x CP^~^ and po = 
(po7 ^^0)7 so that Po is the minimum of V : P ^ M.. Let 
p{t) be any smooth curve in P with p{0) = po^ and let 
p{0) = X e Tp^P. Since po is a critical point of 
dVp^X = {V op)'(O) = 0. Now HesSpQ is, by deflnition, 
the unique symmetric bilinear form on Tp^P such that 



d'V{p{t)) 



Hess^,(X,X) 



(5.30) 



^=0 



for all curves p{t). Since po is a minimum of HesSpQ 
is non-negative, that is, HesSpQ(X, X) > for all X. The 
vector space Tp^P is equipped with an inner product. 



{{cT,Y),{a\Y^)h 



1 



(5.31) 

SO we can uniquely identify Hess^Q with a self- adjoint lin- 
ear map HpQ : Tp^P TpoP such that 



HesSp,{X,X') = {X,np,X'). 



(5.32) 



This is the energy functional of a set of decoupled flelds, 
consisting of a Proca (vector boson) fleld J of mass 



mj epo 



(5.35) 



and {2k — 1) real Klein-Gordon (scalar boson) flelds a^, 
of masses m^. 

In general, the asymptotic flelds of a soliton will have 
all these degrees of freedom non-zero, and the dominant 
force between well-separated solitons will be mediated 
by whichever mode has longest range, that is, lowest 
mass. So the flrst task in predicting long range intersoli- 
ton forces is to compute the spectrum of the self-adjoint 
linear map ^(po,<^>o)- ^ generic choice of V in the 
family we are considering (2.1), it is not possible to com- 
pute even the vacuum (po^^o) explicitly, so the matrix 
^{po,^o)^ and hence its spectrum, is perforce known only 
numerically. There are, however, some interesting cases 
where explicit analytic progress is possible. 



D. The sigma model limit 

In this section we consider the /c-component GL model 
with potential 



V=-A(l-|^ 



2\2 



(5.36) 



in the limit A ^ 00, where 77 is a real-symmetric k x 
k matrix, with zero diagonal, parametrizing a general 
collection of Josephson interactions. In the notation of 
Sec. II this is the case — = Pa = lab = A for all a, b. 
The special case where 77 = 0, A ^ 00 and e ^ 00, 
which reduces to a pure sigma model, was considered in 
Refs. 25 and 26. It is possible to flnd explicit formulae for 
the topological solitons in that case. The case of flnite A 
and e, with 77 = 0, has also been treated previously. ^^^^ 
The fleld equations for the model (5.25) in the sigma 
model limit (in fact, in the case where ^ is valued in any 
compact kahler manifold) were studied in detail, from a 
geometric viewpoint, in Ref. 30. Our focus here is on the 
new phenomena introduced by the Josephson terms rj. 
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Figure 16. (Color online) - The Q — 1 soliton for the t/(3) symmetric model broken by Josephson interactions of the form 



(5.42), with A = 20 and ryo = 1. The quantities |1 



(panel A), and lej^l (panel B), measure the deviation from the 



cr-model. They converge to zero as A is increased. Panel C shows the energy density of the skyrmion.The fourth panel (panel 
D) displays the texture of the field n, which is similar to that of a baby skyrmion. 



In terms of the polar coordinates p, Z, the limit A ^ oo 
amounts to the constraint p = 1, and the potential V 
reduces, in this limit, to 



vm) 



1 Z^r^Z 



(5.37) 



We have included the factor of |Zp in the denominator of 
this expression (which, of course, equals 1 since \Z\ = 1 
by definition) so that the right hand side is manifestly 
a function of the projective equivalence class of Z only, 
not Z per se, that is, V{[cZ]) = V{[Z]) for all c G C\{0}. 
This is convenient when one comes to compute the Hes- 
sian of V. Since rj is real symmetric, it has a unitary ba- 
sis of eigenvectors ei, 62, . . . , e/c, with corresponding real 
eigenvalues Ai > A2 > • • • > A/^. Expanding Z relative to 
this basis 



k 

E 

i=l 



we see that 



V 



k 



(5.38) 



(5.39) 



Hence, the U{k) symmetry of the model, which is pre- 
served by the sigma- model limit, is broken by rj generi- 
cally to [/(l)^. In the case where the spectrum of r] is 
degenerate, the breaking may be partial. For example, 
if Ai = A2 and all other A^ are distinct, the free energy 
remains invariant under /7(2)x[/(l)^~^, where U{2) acts 
in the obvious way on the span of {ei, 62}. 

Clearly, V : CP^"^ R attains its minimum at 
[Z] = [ei], and this minimum is unique if Ai ^ A2. If 
Ai = A2 = • • • = Aj > Aj+i > • • • > A/e, then any Z in 
the span of {ei, . . . , e^} minimizes so the set of min- 
ima of y is a CP^-^ submanifoldof CP^-^ In this case, 
there can be no energy minimizer on M? with Q ^ 0, by 
Derrick's scaling argument, {i.e. solitons are unstable 
against expanding indefinitely) so let us assume, hence- 
forth, that Ai 7^ A2, so that the vacuum of the model. 



[ei], is unique. If the field <!> = tt o ^ : ^ CP^"^ has 
topological charge Q = I then it wraps R^ once around 
some submanifold homologous to CP^ in CP^~^. In or- 
der to minimize the contribution of it should be the 
CP^ on which ^ lies in the span of {61,62}, the sum of 
the two highest eigenspaces of 7^. So we predict that 



(5.40) 



everywhere, where Xi = complex valued func- 

tions on R^. From the pair (X17X2) we can construct a 
5'^-valued field using the usual identification of CP^ with 
5*^, that is 



n = (Xi X2)^ 



Xi 

X2 



(5.41) 



where r = (ti, r2, ra) are the Pauli spin matrices. In this 
way, a Q = 1 energy minimizer can, conjecturally, be 
identified with a degree 1 texture n : R^ ^ 6*^. Since ^ 
is parallel to 61 at |x| = 00, we see that X2(oo) = 0, and 
hence 11(00) = (0,0,1)^. 

We present numerical evidence in favor of this conjec- 
ture in Fig. 16, in the case /c = 3, 



-1 -1 
r] = Vo\ -I -2 
-1 -2 



(5.42) 



770 = 1 and A = 20. It is found that ^ approximately 
satisfies the sigma- model constraint, more precisely, Si = 
maxa,^^2 |1 — l^pl < 0.04. For this choice of r], 



Ai = 2r]o , ei = 



A2 (V^- l)/?o ,62 



A3 -(V^+ 1)7?0 ,63 



1 








' ( 


' v^+l 




+ 2V3 \ 




' ( 






-2V3 \ 
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We expect the Q = 1 energy minimizer to have ^ in 
the span of {61,62} which, since the eigenvectors form 
a unitary frame, is equivalent to satisfying 63^ = 0. 
Again, this turns out to be approximately true: £2 = 
max2;^^2 163^1 < 0.03. We find that both the errors 
ei and S2 become smaller as A increases with 770 held 
fixed. This indicates that the sigma model limit is well 
founded and should be a reliable approximation for A 
large but finite. Qualitatively, in this special case of the 
3-component model, the n field we find numerically is 
similar to the field of a so-called baby-skyrmion.^^ 

If we place two Q = 1 energy minimizers a long dis- 
tance apart and allow the system to relax, do they repel 
one another and escape to infinity, or do they attract 
one another and coalesce into a Q = 2 bound state? To 
predict this, we need to compute the spectrum of the 
Hessian of V about <l>o = [61], as described in Sec. VC. 
In this case, p is frozen by the constraint, so P = CP^~^. 
It is useful to identify the tangent space Tfg^jCP^""*^ with 
the {k — l)-dimensional complex vector space 



= {Y eC^ : e\Y = 0}. 



(5.43) 



Then the natural metric on Tp^P (5.31) reduces to 

{Y,Y')v = Re{Y^Y') (5.44) 

the restriction of the Euchdean metric on C'' to V. To 
compute the Hessian of V about [ei], we consider a curve 
Z{t) in C'^ with Z{0) = d and Z{0) = F G V. Then 



Hess[e,](y,r) 



d- 



(5.45) 



where we have used the fact that Xilk — r] is self adjoint 
with 6i in its kernel, so that e\{XIk — r])Y = 0. Hence, the 
associated self-adjoint linear map H[e^] : V ^ V is the 
restriction to V of A il^ —77. It follows that the eigenvalues 



of 1-L[e^] are Ai — A^, i = 2, 



, /c, each of multiplicity 



2, and that the corresponding eigenspaces are two real- 
dimensional, spanned by {6^, ici}^ z = 2, . . . , /c. So there 
are 2k — 2 real scalar bosons in this model, occurring in 
pairs, having mass 



rrii 



VAi-A,. 



This should be compared with the mass of the supercur- 
rent field, i.e. the inverse London penetration length. 



mj 



(5.47) 



Numerics suggest that the supercurrent of a Q = 1 
energy minimizer is, at large similar to that of a 
vortex, while the lightest (complex) Klein-Gordon mode 
X2 is similar to the asymptotic field of a baby-skyrmion. 
Hence, we expect J to mediate a repulsive force of range 



1/6 and X2 to mediate a short-range scalar dipole-dipole 
force. The range of this force is l/^/Xl^^J^. The lat- 
ter force is attractive provided the two solitons are ap- 
propriately aligned; see the discussion of baby-Skyrme 
models^^ for a detailed analysis. The dipole like interac- 
tion is also natural from the viewpoint of the fractional- 
vortex picture of skyrmions (see discussion in Sec. IV and 
in Refs. 18). Hence, we predict that a pair of Q = 1 soli- 
tons, in the model which we consider in this subsection, 
always repel (for all relative orientations) if 6^ < Ai — A2, 
so higher Q bound states cannot form. On the other 
hand, if > Ai — A2, well-separated solitons have an 
attractive channel, and we predict that they can coa- 
lesce into higher Q bound states. Numerical evidence 
of this predicted dichotomy in the three component case 
is presented in Fig. 17 and direct numerical evidence of 
dipolar interaction of two Q = 1 skyrmions is presented 
in Fig. 18. 



E. A symmetric case with BTRS 

In this section we consider the GL^^^ model with three 
identical active bands, coupled through identical Joseph- 
son terms. The potential is 

!/(*) = ^ E(l - + i*^^*' (5-48) 

0=1 

where N denotes the symmetric coupling matrix 

TV = I 1 1 1 . (5.49) 



Note that, in contrast to section VD, 7^ denotes a real 
parameter here, not a matrix. In terms of the notation 
of section II, this is the special case ai = 0^2 = 0^3 = — ^, 

Pi = p2 = p3 = lab and r]i2 r]i3 V23 -V- 
The vacuum manifold for this potential is a disjoint union 
of two circles, the gauge orbits of 




^ = pqVq and ^ = poVi 



where 



Po = a/3 



6r] 
T' 



(5.50) 



(5.51) 



(5.46) and (with ^ = e^^'/'') 




^ (5.52) 

are simultaneous unit eigenvectors of the symmetric cou- 
pling matrix N and the permutation matrix P, 



(5.53) 
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Figure 17. (Color online) - The results of energy minimization with charge Q = 2 for A = 20, and e = 1.0 and varying 770. First 
row shows the energy density while the second and third row displays the corresponding texture field. Panels (A, 770 = 0.1), 
(B, rjo = 0.2)), (C, ?7o = 0.3)), are for 770 < tjq where interaction between skyrmions is attractive. Two Q — 1 skyrmions 
coalesce into either one Q = 2 skyrmion (A, B and C). Configuration displayed on panel C resemble bound state of two Q = 1 
skyrmions. Panel (D, 770 = 0.8)) has r]o > r]Q, then in the repulsive channel. Here the two Q = 1 skyrmions are repelling 
each other. So the snapshot on panel D shows a late but unconverged iteration (i.e. it represents a fairly converged pair of 
individual skyrmions which are, however, still drifting apart). 
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where 



Figure 18. (Color online)- Interaction energy of two single 
quantum solitons. The the GL parameters are aa — —20, 
Pa = 20, 7ab = 20, Tjah = — 1 and e = 1. Thus, the potential 
part of the free energy density can be written as [/ = A(l — 

iV^ll^ - |'02|^ - iV^sl^)^ - ryabl^allV^bl COs((/?a " ^h) with A = 

10, i.e. the Hamiltonian features SU{3) symmetry broken by 
Josephson interaction term. 



Note that 



-vo, Nvi 



Nv2 = 2v2, (5.54) 



(5.55) 



We shall, without loss of generality, choose the vacuum 
Po'^o (rather than po'^i)- Since [^] 7^ [^] for this vacuum, 
the model has broken time-reversal symmetry. 

There are axially symmetric vortex solutions which in- 
terpolate between (0,0,0) at r = and the above vacuum 
at r = 00. To construct them, one only needs to solve a 
single component GL model: 




{pI-\<p\ 



2\2 



Given a vortex solution (^, A) of (5.56), 



(5.56) 



(5.57) 



is a vortex solution of the symmetric GL^^^ model (5.48). 
The numerical results of Sec. II strongly suggest that 
(5.48) also supports skyrmion solutions, at least for Q 
and 1/e sufficiently large. 

Once again, we wish to compute the spectrum for the 
Hessian of V about the vacuum (po, bo])- The potential 
is, in polar coordinates (5.12), 



V{p, Z) = ^ ^(1 - p'\ZXf + Ip'Z^NZ (5.58) 

a=l 

= Y - + -/um + Imz]) (5.59) 



urn) 



' ' a=l 



(5.60) 
(5.61) 



We have included the factors of |Zp in the denominators 
of these expressions (which, of course, equals 1 by defini- 
tion) so that the right hand sides are manifestly functions 
[Z] only. Recall that Hess is a symmetric bilinear form 
on the tangent space to (0, 00) x CP^ at the vacuum 
(po, bo])- In general, there is no reason why this bilinear 
form should not couple the direction tangent to (0, 00) 
with directions tangent to CP^ . We shall see that in this 
case permutation symmetry prevents such coupling. 

First, we note that [vq] is a fixed point of the permu- 
tation map 



V : CP^ ^ CP^ 
and that dPi 



[Z]^[PZ], 



(5.62) 



: T[^q]CP^ T[^q]CP^ has maximal rank, 
so it follows that [vq] is a critical point of any function 
CP^ M invariant under P. In particular. 



0. 



(5.63) 



Consider now a two-parameter variation p{s,t) = 
{p{s), [Z{t)]) through po = (po, N]) in P = (0, ^) xCP^, 
with 5sp(0,0) = ((j,0) and dtp{Q,0) = (0,r). Then 



HesSp„((cT,0),(0,r)) 



d^V{p{s,t)) 



dsdt 



s=t=0 



pladU[^^]Y + r]poadU[y„]Y 

(5.64) 



by (5.63). Hence 

Hess = (A + 2ri)dp- 



(5.65) 

where U. il : Ty^^\CP^ x Tf^ojCP^ R are the Hessians 

of the functions [/, U respectively. It follows that one of 
the real scalar bosons ai in (5.34) is just a (the lineariza- 
tion of p about po) and that this has mass 



Vflp = VA + 27/. 



(5.66) 



It remains to compute H and R. For this purpose, we 
identify the tangent space T[^q]CP^ with the two dimen- 
sional complex vector space 



V = {r G : vIy = 0} 



(5.67) 



which is spanned by {v\^V2\^ and give V the induced 
Euclidean metric 

{X,F)v= ^(Xtr + ytx) = (5.68) 
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where (•,-)fs' denotes the Fubini-Study metric, used to 
compute IdFp in equation (5.25). 

In fact, we aheady know since this is a special case 
of the general Josephson coupling matrix considered in 
section VD: 



H{X,Y)=2{X, (7V + l)rK 



(5.69) 



It is convenient to expand X, Y relative to the unitary 
(for basis Vi^V2^ which are eigenvectors of N. 

Namely, if 



X = (xi + iX2)Vi + (X3 + ix/^)v2, 
Y = {yi^ iy2)vi + (^3 + m)v2 



then 



H{X,Y) =6x^ 



/O 



10 

V 1 



y- 



(5.70) 



(5.71) 



Note that this is a hermitian bilinear form on V, that is 
H{iX,iY) = H{X,Y). 

Turning to one should not expect it to be hermitian, 

2 

because U contains terms like ZfZ^. Consider a two- 
parameter variation Zs^t in C with Zo,o = '^0 and 
dsZs,t\o,o = X G V, dtZs^t\o,o = V. By definition. 



H{X,Y) = 



d^U{Zs,t) 



dsdt 



(5.72) 



=t=0 



Using the explicit formula (5.60) for U{Z)^ we find that 



H{X,Y)=2j2i^aZa^ZaXa){YaZa^ZaYa) (5.73) 



a=l 



where Z = vq. Note this is not Hermitian because, for ex- 
ample H{ivi^iv2) = —H{vi,V2). Again, we can express 
this as a 4 X 4 real matrix, by expanding X, Y relative to 
vi^V2. One finds that 



H{X,Y) = -x' 



/ 1 1 
10-1 
10 10 

Vo -1 1 



y- 



(5.74) 



Substituting (5.74) and (5.71) into (5.65), then (5.65) 
into (5.29), we obtain 



(iidjr+e^p^iijr) 
1 



+ ^i\da\' + (A + 2n)a') + -pi (|dt/|^ + My) 



where the mass matrix is 



M = 



/ 1 1 
10-1 
10 10 

Vo -1 1 



-3r] 



/O 

10 

Vo 1 



(5.75) 



. (5.76) 



The squared masses of the bosons tangent to CP^ are 
the eigenvalues of this matrix, namely 




(5.77) 



each of multiplicity two. These should be compared with 
the mass of the J vector boson and p scalar boson 



m 7 



2 2 



(5.78) 



To extract information about intersoliton forces, note 
that the embedded vortex (5.57) excites only the (repul- 
sive) J mode and the (attractive) p mode, so one pre- 
dicts the usual behaviour {i.e. for the example consid- 
ered here where there is degeneracy in couplings between 
components, at long range vortices repel if rrip > mj, 
and attract if rup < mj). Note that in the case when 
the components have different prefactors in V, there 
are also type- 1.5 regimes with non-monotonic intervor- 
tex (long-range attractive, short-range repulsive) inter- 
vortex forces.^ Skyrmions, on the other hand, should in 
all cases excite all 6 modes, with a monopole source for p 
and dipole (or higher) sources for the 4 (mixed) Y modes. 
So an interesting regime would be m_ < mj < nip since 
then inter vortex forces should be long-range repulsive, 
while inter-skyrmion forces should have an attractive 
channel for a certain relative orientations of skyrmions. 



VI. CONCLUSIONS 

We discussed a new kind of topological soliton which 
we term chiral GL^^^ skyrmions. These solitons oc- 
cur in three-component superconductors when time re- 
versal symmetry is spontaneously broken. In contrast 
to vortices, these skyrmions are characterized by a 
(DP^ topological charge. These skyrmions have a defi- 
nite chirality associated with them: i. e. the order of the 
constituent fractional vortices matters, different orders 
giving inequivalent solutions. We described two situa- 
tions 

• A type-II BTRS superconductor can form a vor- 
tex lattice as a ground state in applied magnetic 
field. However in contrast to usual vortex states, 
all the regimes investigated by us possessed other 
fiux-carrying topological defects of a higher en- 
ergy: metastable GL^^^ skyrmions characterized 
by a (DP^ topological charge. The system thus 
can form infinitely many complex metastable states 
in external fields where vortices coexist with the 
GL^^^ skyrmions solitons. Thermal, or magnetic 
field quench can force the system to fall into one of 
these states. 

• BTRS three-band superconductors in principle can 
have also a different regime where in external field 
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(DP^ solitons are energetically cheaper than vor- 
tices. In that case the system cannot form vor- 
tices since they are unstable against decay into 
skyrmions. Such regimes occur for example when 
the free energy has bi-quadratic interaction terms 
of the form 7a6|V^aPlV^6p. 

In the regimes where chiral GL^^^ skyrmions are 
metastable they can spontaneously form from 'collisions' 
of vortices, where intervortex interaction energy can be 
larger than energy of potential barrier of forming a 
skyrmion. We investigated several hundred regimes and 
found that skyrmions typically easily form in the energy 
minimization process where a system is relaxed from var- 
ious higher energy states (such as dense groups of ordi- 
nary vortices). Our study indicates that the "capture 
basin" of these solutions can in certain cases be very 
large. We find that these defects very easily form dur- 
ing a rapid expansion of a vortex lattice (which should 
occur when magnetic field is rapidly lowered, or if a sys- 
tem is quenched through Hc2)- Formation of solitons in 
this process can signal a state with Broken Time Re- 
versal Symmetry. Also the potential barriers between 
Skyrmions and vortices or between different skyrmionic 
states can be overcome due to thermal fiuctuations. 

As shown in Fig. 1, these skyrmions have very par- 
ticular magnetic signature and thus, under certain con- 
ditions, may be observed in high-resolution scanning 
SQUID, Hall, or magnetic force microscopy measure- 



ments. A tendency for vortex pair formation, yield- 
ing magnetic profile similar to that shown on Fig. 5 
was observed in Ba(Fei_xCox)2As2,^^ as well as vor- 
tex clustering in BaFe2-xNixAs2.^^ These materials have 
strong pinning which can naturally produce disordered 
vortex states, although the possibility of "type- 1.5" sce- 
nario for these vortex inhomogeneities was also voiced 
in Ref. 35 . (Note that in three band (or higher num- 
ber of bands) superconductors with frustrated Joseph- 
son coupling, type- 1.5 regimes are easily obtainable even 
if Josephson coupling is very strong.^) The vortex pairs 
observed in Ref. 34 can be discriminated from Q = 2 
solitons by quenching the system in a stronger magnetic 
field and observing whether or not it forms vortex tri- 
angles, squares, pentagons, such as shown on e.g. Fig. 1 
which correspond to fiux profile of higher- Q solitons. Be- 
sides multiband superconductors, another class of sys- 
tems which can support chiral GL^^^ skyrmions is a 
Josephson coupled sandwich of an s± and 5-wave super- 
conductor. 
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Appendix A: Fractional vortices have linearly 
divergent energy in the presence of Josephson 
coupling 

Here we discuss fractional flux vortices in three band 
systems. Consider the case of one fractional vortex in 
which winds through 27r and neither ?/^2 nor ?/^3 winds. 
We assume that the configuration is spatially localized 
around r = 0, so that on any annulus Q = {ro < r < ri}, 
with ro sufficiently large, the densities {ipal are close to 
their ground state values {i.e. we assume the London 
limit). It follows from expression (2.1) that the total 
free energy of any configuration satisfies the lower bound 



F - Fo. > R 



GS ^ SG 7 



(A.l) 

with FsG''=^J^ab 1 \\/(Pab\'^ ^l,rnl^{^ - cos (fab) 
a<b ^ 

where Vab = |V^aPlV^6p/p^, = "^Vabp^ /\^a\\'4^b\, and 
= IV^aP- ^Gs denotes the energy of the vortex- 
less ground state. In the London limit, the field densities 
assume their ground state values, so Uab and rriab are 
constants. In this limit, F^g simplifies to a sum of sine- 
Gordon energies (hence the subscript SG). Note that 
|V(/:'a6p ^ r~'^{d(pab/d6)'^ ^ with r and ^, the polar coor- 
dinates around the vortex center. Hence, 



FsG > ^a6 / \ 
a<b I 

a<b I 



r2 V do 

1 d(pab 

r dO 



2 • 2 Yab 



(A.2) 



rriab sm - 



^ab 



2mab difab . ^ab 



> 



(A.3) 



a<b "^^0 

8(^12^12 + rnisiyis){ri - ro) 



(A.5) 



where we have used the boundary conditions that (pi2 
and (fis wind once, while (/?23 does not wind. So Fgc, 
and hence the total free energy F — Fqs, grows (at least) 
linearly with the system size, ri. 

Note that our lower bound on Fgo cannot be attained, 
because for this to happen, one would need (fab to satisfy 



1 d^ab 

r dO 



rriab sm - 



^ab 



(A.6) 



and no solutions to this PDE with the correct boundary 
behaviour ((/:?i2(r, 27r) — (pi2{r,0) = 27r for all r) exist. 



Appendix B: Finite element energy minimization 

The chiral skyrmions are either global or local minima 
of the Ginzburg-Landau energy (2.1). In the later case. 



this means that a good enough initial guess is necessary. 
In both cases, the functional minimization of (2.1), from 
an appropriate initial guess carrying several flux quanta, 
should lead to a chiral skyrmion (if it exists as a sta- 
ble solution). We consider the two-dimensional problem 
(2.1) defined on the bounded domain Q C with 
its boundary. In practice we choose ^7 to be a disk. Ac- 
tually, the particular shape of the domain is not impor- 
tant. Indeed it is much larger than the typical size of 
solitons. Moreover, neither solitons nor initial guess co- 
incide with some grid symmetry. For example, skyrmions 
are never placed at the center of the domain (the vizual- 
ization scheme re-centers the window around the soli- 
ton). This is an addtional argument that skyrmions are 
not boundary artefacts. One some occasions, we dou- 
bled checked on square domains that our solutions are 
unaffected by boundaries. 

The problem is supplemented by the boundary con- 
dition n ■ D^a = with n the normal vector to dVt. 
Physically this condition implies there is no current flow- 
ing through the boundary. Since this boundary condition 
is gauge invariant, additional constraint can be chosen on 
the boundary to fix the gauge. Our choice is to impose 
the radial gauge on the boundary • A = (note that 
with our choice of domain, this is equivalent to n- A = 0). 
With this choice, (most of) the gauge degrees of freedom 
are eliminated and the 'no current flow' condition sepa- 
rates in two parts 



and 



n- A = 0. 



(B.l) 



Note that these boundary conditions allow a topological 
defect to escape from the domain, since there is no pres- 
sure of an external applied field. Because they are topo- 
logical defects, vortices (and skyrmions) cannot unwind. 
However, they can be 'absorbed' through the boundary 
in order to further minimize the energy. To prevent this, 
the numerical grid is chosen to be large enough so that 
the attractive interaction with the boundaries is negligi- 
ble. The size of the domain is then much larger than the 
typical interaction length scales. Thus in this method one 
has to use large numerical grids, which is computation- 
ally demanding. The advantage is that it is guaranteed 
that obtained solutions are not boundary pressure arti- 
facts. 

The variational problem is defined for numerical com- 
putation using a finite element formulation provided by 
the Freefem++ library. Discretization within finite el- 
ement formulation is done via a (homogeneous) trian- 
gulation over 1^, based on Delaunay-Voronoi algorithm. 
Functions are decomposed on a continuous piecewise 
quadratic basis on each triangle. The accuracy of such 
method is controlled through the number of triangles, 
(we typically used 3 ~ 6 x 10^), the order of expansion 
of the basis on each triangle (2nd order polynomial basis 
on each triangle), and also the order of the quadrature 
formula for the integral on the triangles. 

Once the problem is mathematically well defined, a 
numerical optimization algorithm is used to solve the 
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Figure 19. (Color online) - A Q = 3 quanta soliton in a system with three identical passive bands as in Fig. 3, except 
that there are no density-density interactions 7ab = and e = 0.3. Since the three bands are identical, the soliton makes a 
homogeneous ringlike configuration. Displayed quantities are the same as in rest of the paper. 

variational nonlinear problem {i.e. to find the minima of 
F). We used here a nonlinear conjugate gradient method. 
The algorithm is iterated until relative variation of the 
norm of the gradient of the functional F with respect to 
all degrees of freedom is less than 10~^. 



Initial guess for obtaining metastable configurations 

As discussed in the paper, N quanta chiral skyrmions 
can be more energetically expensive than N ordinary 
(type-II) vortices. In that case the initial guess should 
be within the attractive basin of the chiral skyrmions. 
Otherwise the configuration converges to ordinary type- 
II vortices which have the same total phase winding but 
cost less energy. The initial field configuration carrying 

fiux quanta is prepared by using an ansatz which im- 
poses phase windings around spatially separated A^ vor- 



tex cores in each condensates. 



\i'a\=Ua 11^11(1 + tanh y)-U 



k=l 



(B.2) 



where a = 1,2,3 and Ua is the ground state value of each 
condensate density. The parameters parametrize the 
core size while 

N 

' y-Vk 



^a{x,y) = ^tan ^ 



X — Xi 



ni{x,y) = ^{x-xlY + {y-ylY. (B.3) 

(x^, determines the position of the core of /c-th vortex 
of the a-condensate.The functions A^^ are used to seed a 
domain wall. As an initial guess we generally choose 
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Figure 20. (Color online) - A Q = 4 quanta soliton in a system with two identical passive bands as in Fig. 5 coupled to a 
third active band with disparity in the ground state densities (as^fis) = (—1-5, 1). Josephson coupling constants are 7723 = —3 
and ?7i2 — 7713 = 1. e = 0.2 and 'jab — 0. 



Ao = —A3 = A, with A defined as 



ro)-l) , 



(B.4) 



where H{t — tq) is a Heaviside function. Thus in the 
initial guess the domain wall has infinitesimal thickness. 
It takes only a few steps from this initial guess to relax 
to a true domain wall during the simulations. Conse- 
quently, it is entirely sufficient to use Heaviside functions 
for the initial guesses for domain walls. The starting con- 
figuration of the vector potential is determined by solv- 
ing Ampere's law equation of (2.2) on the background 
of the superconducting condensates specified by (B.2)- 
(B.4). Being a linear equation in A, this is an easy op- 
eration. 

Once the initial configuration defined, all degrees of 
freedom are relaxed simultaneously, within the 'no cur- 
rent flow' boundary conditions discussed previously, to 
obtain highly accurate solutions of the Ginzburg-Landau 
equations. In a strongly type-II system when the ini- 



tial guess was either (a) vortices placed on a closed 
domain wall or (b) closed domain wall surrounding a 
densely packed group of vortices, the system almost 
always formed chiral skyrmions. We used also initial 
guesses (c) without any domain walls (A = 0). In that 
case we observed chiral skyrmion formation, if in the ini- 
tial states vortices were densely packed. This again indi- 
cates that the chiral skyrmions in the three component 
GL model represent (local) minima with wide capture 
basin in the free energy landscape. 



Appendix C: Additional Material 

In this appendix we show few additional solutions 
Fig. 19, Fig. 20 and Fig. 21 for chiral skyrmions. Pa- 
rameters sets, or number of flux quanta used here are 
different from the ones considered in the main body of 
the paper. 
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Figure 21. (Color online) - A Q — 7 quanta soliton in a system with two identical passive bands as in Fig. 5 coupled to a 
third active band with disparity in the ground state densities (asj/^s) = ( — 1, 1). Josephson coupling constants are 7723 = —3 
and ?7i2 — 7713 = 1. e = 0.3 and — 0. 



